function [lb,ub] = fsel_pi(X,Y,Xc,set0,kmax)
% Forward selection for inference for profit differential functional
% Last edit: 3/28/17 CBH
% [lb,ub] = fsel_pi(X,Y,Xc,set0,K)
% X - design matrix
% Y - outcome
% Xc - matrix of evaluation points for profit functional
% set0 - set of initial regressors
% kmax - number of forward selection steps to take

p = size(X,2);
varind = (1:p)';

supp_lb = set0;
supp_ub = set0;

lb = zeros(kmax,1);
ub = zeros(kmax,1);

for jj = 1:kmax
    varind_lb = setdiff(varind,supp_lb);
    varind_ub = setdiff(varind,supp_ub);
    
    ci_lb = zeros(numel(varind_lb),1);
    ci_ub = zeros(numel(varind_ub),1);
    
    for kk = 1:numel(varind_lb)

        use_lb = [supp_lb;varind_lb(kk)];
        zlb = X(:,use_lb);
        b_lb = zlb\Y;
        e_lb = Y - zlb*b_lb;
        BLB = zeros(p,1);
        BLB(use_lb) = b_lb;

        use_ub = [supp_ub;varind_ub(kk)];
        zub = X(:,use_ub);
        b_ub = zub\Y;
        e_ub = Y - zub*b_ub;
        BUB = zeros(p,1);
        BUB(use_ub) = b_ub;
        
        ci_lb(kk,1) = TUSIM_prof(BLB,Xc) - 1.96*TUSIM_profse(BLB,Xc,1e-6,X,use_lb,e_lb);
        ci_ub(kk,1) = TUSIM_prof(BUB,Xc) + 1.96*TUSIM_profse(BUB,Xc,1e-6,X,use_ub,e_ub);
        

    end
    imin = find(ci_lb == min(ci_lb));
    imax = find(ci_ub == max(ci_ub));
    imin = imin(1);
    imax = imax(1);
    
    supp_lb = [supp_lb ; varind_lb(imin)]; %#ok<*AGROW>
    supp_ub = [supp_ub ; varind_ub(imax)];
        
    lb(jj,1) = min(ci_lb);
    ub(jj,1) = max(ci_ub);
end
